Eel silvering model summary

Author

Viktor Thunell

Published

January 22, 2026

Load libraries

Show code
# Load libraries, install if needed
pkgs <- c("tidyverse", "tidylog", "devtools","viridis","sf","patchwork","coda","boot","tidybayes","bayesplot", "nimble", "here","reshape2")

if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){
    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  }

invisible(lapply(pkgs, library, character.only = T))

options(ggplot2.continuous.colour = "viridis")
theme_set(theme_minimal()) 

# Set path
home <- here::here()

Read data

Show code
## Read data
df.eel <- readRDS(file = "~/gitProjects/2024_R_DIASPARA/data/indhdb_3.RData") %>%
  st_drop_geometry()

Length at Silvering

Define data

Show code
dfp <- df.eel %>%
  filter(fi_lfs_code == "S",
         fi_year > 1997,
         # norway has three silver eels in total, so I remove this
         !ecoreg == "norwegian") %>% 
  rename(length = lengthmm) %>%
  mutate(sex = as.integer(if_else(sex == "f", 0, 1)))
filter: removed 724,965 rows (83%), 147,494 rows remaining
rename: renamed one variable (length)
mutate: converted 'sex' from character to integer (0 new NA)
Show code
data.silv <- dfp %>%
  mutate(mb = as.integer(factor(main_bas)),
         er = as.integer(factor(ecoreg)),
         year = as.integer(factor(fi_year, levels = sort(unique(fi_year)), labels = seq_along(sort(unique(fi_year))))),
         temp.sc = (dfp$tmp_dc_uyr - mean(dfp$tmp_dc_uyr))/sd(dfp$tmp_dc_uyr),
         age.sc = (dfp$ageyear - mean(dfp$ageyear))/sd(dfp$ageyear)) 
mutate: new variable 'mb' (integer) with 154 unique values and 0% NA
        new variable 'er' (integer) with 9 unique values and 0% NA
        new variable 'year' (integer) with 28 unique values and 0% NA
        new variable 'temp.sc' (double) with 78 unique values and 0% NA
        new variable 'age.sc' (double) with one unique value and 100% NA
Show code
rm(dfp)

Plot silvering data

Show code
# length distributions
data.silv %>%
  mutate(sex = if_else(sex == 0, "f", "m")) %>%
  ggplot() +
  geom_density(aes(x = length, fill = sex), alpha = 0.5) +
  facet_wrap(~ecoreg) +
  theme_minimal()
mutate: converted 'sex' from integer to character (0 new NA)

Show code
data.silv %>%
  mutate(sex = if_else(sex == 0, "f", "m")) %>%
  ggplot() +
  geom_density(aes(x = log(length), fill = sex), alpha = 0.3) +
  labs(title = "log(length)") +
  facet_wrap(~ecoreg) +
  theme_minimal() 
mutate: converted 'sex' from integer to character (0 new NA)

Show code
# length by age
data.silv %>%
  filter(!is.na(ageyear)) %>%
  mutate(sex = if_else(sex == 0, "f", "m")) %>%
  ggplot() +
  geom_point(aes(ageyear, length, color = sex), size = 0.1, alpha = 0.8) +
  facet_wrap(~ecoreg, scales = "free_y")  +
  theme_minimal()
filter: removed 131,289 rows (89%), 16,205 rows remaining
mutate: converted 'sex' from integer to character (0 new NA)

Show code
# possible temperature effects on silvering length
data.silv %>%
  mutate(sex = if_else(sex == 0, "f", "m")) %>%
  ggplot(aes(tmp_dc_uyr/10, length)) +
  geom_point(show.legend = FALSE, alpha = 0.5) +
  geom_smooth(method = "lm", se = TRUE) +
  facet_wrap(~factor(sex)) +
  theme_minimal() +
  labs(x = "Temp. [C]")
mutate: converted 'sex' from integer to character (0 new NA)
`geom_smooth()` using formula = 'y ~ x'

Source models / Load samples

Show code
# Load samples
# readRDS(paste0(home,"/data/silv.samples_1214.RData"))$WAIC
# readRDS(paste0(home,"/data/silv.samples_1216.RData"))$WAIC
silv.samples <- readRDS(paste0(home,"/models_eel/samples/silv.samples_2026-01-08.RData"))$samples # v9

# # Or source model code and mcmc (load libs and data)
# t <- Sys.time()
# source(file = paste0(home,"/Eel models/2-2_model_eel_silver_v5.R"))
# Sys.time() - t # approx xxx hours with xx million hmc samples

Traces, ess, and ‘potential scale reduction factor’

Variables “s.sig”,“Ustar”,“U”,“mu.a.gl”,“sd.a.gl”,“alpha”,“bs”,“bt”)

Show code
gc()
            used   (Mb) gc trigger   (Mb) limit (Mb)  max used   (Mb)
Ncells   5102104  272.5   11934382  637.4         NA  11934382  637.4
Vcells 665703217 5079.0  960990790 7331.8     102400 665933856 5080.7
Show code
## bs s.sig bt ba
node.sub <- grep("^bs|s.sig|bt|ba", colnames(silv.samples[[1]]), value = TRUE)
silv.samples[, node.sub[sample(1:length(node.sub))], drop = FALSE] %>%
mcmc_trace()

Show code
silv.samples[, node.sub[], drop = FALSE] %>%
gelman.diag()
Potential scale reduction factors:

      Point est. Upper C.I.
bs[1]      1.000      1.000
bs[2]      1.002      1.003
bs[3]      0.999      0.999
bs[4]      1.000      1.000
bs[5]      1.000      1.000
bs[6]      1.000      1.004
bs[7]      1.001      1.005
bs[8]      1.000      1.002
bs[9]      1.038      1.157
bt         1.014      1.064
s.sig      1.000      1.001

Multivariate psrf

1.03
Show code
print("effective sample size")
[1] "effective sample size"
Show code
silv.samples[, node.sub[], drop = FALSE] %>%
effectiveSize()
    bs[1]     bs[2]     bs[3]     bs[4]     bs[5]     bs[6]     bs[7]     bs[8] 
2405.2796 3000.0000 3000.0000 2474.9394 2309.7995 2364.5004 2710.7620 2397.2730 
    bs[9]        bt     s.sig 
 354.9447  501.0123 2404.8789 
Show code
## alpha
node.sub <- grep("^alpha", colnames(silv.samples[[1]]), value = TRUE)
pl <- silv.samples[, node.sub[], drop = FALSE] %>%
gelman.diag()
pl$psrf[which(pl$psrf[, 1] > 1.01),]
             Point est. Upper C.I.
alpha[4, 1]    1.035842   1.157531
alpha[9, 1]    1.017975   1.084183
alpha[4, 2]    1.036431   1.160147
alpha[9, 2]    1.021571   1.098294
alpha[1, 3]    1.014910   1.066775
alpha[3, 3]    1.011207   1.053803
alpha[4, 3]    1.034445   1.152616
alpha[9, 3]    1.025539   1.114275
alpha[1, 4]    1.013031   1.058254
alpha[3, 4]    1.013710   1.065541
alpha[4, 4]    1.038463   1.167569
alpha[9, 4]    1.023603   1.103274
alpha[1, 5]    1.013582   1.062121
alpha[3, 5]    1.012638   1.058960
alpha[4, 5]    1.038580   1.166649
alpha[9, 5]    1.024120   1.081205
alpha[1, 6]    1.016550   1.074840
alpha[3, 6]    1.014017   1.068200
alpha[4, 6]    1.045089   1.189158
alpha[9, 6]    1.030717   1.120727
alpha[1, 7]    1.019309   1.081682
alpha[2, 7]    1.011234   1.055826
alpha[3, 7]    1.010594   1.052865
alpha[4, 7]    1.043874   1.185690
alpha[9, 7]    1.025584   1.108955
alpha[1, 8]    1.013931   1.062618
alpha[4, 8]    1.023791   1.102031
alpha[9, 8]    1.026379   1.120304
alpha[1, 9]    1.010812   1.049027
alpha[4, 9]    1.016677   1.075781
alpha[9, 9]    1.026811   1.121061
alpha[1, 10]   1.012026   1.044396
alpha[9, 10]   1.024643   1.101231
alpha[1, 11]   1.015123   1.054735
alpha[9, 11]   1.026405   1.100825
alpha[4, 12]   1.011084   1.030381
alpha[9, 12]   1.037431   1.156163
alpha[9, 13]   1.034936   1.148972
alpha[9, 14]   1.041291   1.178526
alpha[7, 15]   1.015286   1.073806
alpha[9, 15]   1.023161   1.106288
alpha[9, 16]   1.024258   1.103777
alpha[9, 17]   1.035194   1.155509
alpha[9, 18]   1.036126   1.158864
alpha[9, 19]   1.038443   1.166945
alpha[9, 20]   1.025539   1.108965
alpha[9, 21]   1.025173   1.112519
alpha[4, 22]   1.013077   1.064108
alpha[9, 22]   1.048991   1.205513
alpha[4, 23]   1.012483   1.055147
alpha[9, 23]   1.044875   1.187352
alpha[5, 24]   1.010856   1.048908
alpha[9, 24]   1.042571   1.170095
alpha[3, 25]   1.011459   1.053094
alpha[9, 25]   1.023989   1.080862
alpha[9, 26]   1.019970   1.065973
alpha[9, 27]   1.017908   1.079335
alpha[9, 28]   1.033370   1.096555
Show code
silv.samples[, node.sub[sample(1:length(node.sub), 50)], drop = FALSE] %>%
mcmc_trace()

Show code
pl <- silv.samples[, node.sub[], drop = FALSE] %>%
effectiveSize()
print("effective sample size")
[1] "effective sample size"
Show code
pl[which(pl < 500)]
 alpha[1, 3]  alpha[1, 6]  alpha[9, 6]  alpha[9, 7]  alpha[9, 8]  alpha[9, 9] 
    481.4697     493.2682     419.1584     471.3161     477.5335     494.2195 
alpha[9, 10] alpha[9, 11] alpha[9, 12] alpha[9, 13] alpha[9, 16] alpha[9, 18] 
    460.8167     396.6974     427.8152     495.6315     404.4980     423.3880 
alpha[9, 19] alpha[9, 20] alpha[9, 21] alpha[9, 22] alpha[9, 23] alpha[9, 24] 
    436.9864     409.3936     397.4059     402.7409     426.5547     437.3401 
alpha[9, 25] alpha[9, 26] 
    440.0765     439.9288 
Show code
node.sub <- grep("^mu.a.gl", colnames(silv.samples[[1]]), value = TRUE)
silv.samples[, node.sub[], drop = FALSE] %>%
mcmc_trace()

Show code
silv.samples[, node.sub[], drop = FALSE] %>%
gelman.diag()
Potential scale reduction factors:

           Point est. Upper C.I.
mu.a.gl[1]      1.008       1.04
mu.a.gl[2]      1.007       1.02
mu.a.gl[3]      1.010       1.05
mu.a.gl[4]      1.031       1.14
mu.a.gl[5]      1.002       1.01
mu.a.gl[6]      1.001       1.00
mu.a.gl[7]      1.001       1.01
mu.a.gl[8]      0.999       1.00
mu.a.gl[9]      1.014       1.06

Multivariate psrf

1.04
Show code
print("effective sample size")
[1] "effective sample size"
Show code
silv.samples[, node.sub[], drop = FALSE] %>%
effectiveSize()
mu.a.gl[1] mu.a.gl[2] mu.a.gl[3] mu.a.gl[4] mu.a.gl[5] mu.a.gl[6] mu.a.gl[7] 
  593.1940  1024.4422   659.9732   683.5194  2541.6480   847.3592   814.3498 
mu.a.gl[8] mu.a.gl[9] 
 2578.9742   943.1168 
Show code
node.sub <- grep("^sd.a.gl", colnames(silv.samples[[1]]), value = TRUE)
silv.samples[, node.sub[], drop = FALSE] %>%
mcmc_trace()

Show code
silv.samples[, node.sub[], drop = FALSE] %>%
gelman.diag()
Potential scale reduction factors:

           Point est. Upper C.I.
sd.a.gl[1]       1.00       1.00
sd.a.gl[2]       1.00       1.02
sd.a.gl[3]       1.00       1.00
sd.a.gl[4]       1.00       1.02
sd.a.gl[5]       1.00       1.01
sd.a.gl[6]       1.00       1.01
sd.a.gl[7]       1.00       1.00
sd.a.gl[8]       1.00       1.01
sd.a.gl[9]       1.02       1.06

Multivariate psrf

1.02
Show code
print("effective sample size")
[1] "effective sample size"
Show code
silv.samples[, node.sub[], drop = FALSE] %>%
effectiveSize()
sd.a.gl[1] sd.a.gl[2] sd.a.gl[3] sd.a.gl[4] sd.a.gl[5] sd.a.gl[6] sd.a.gl[7] 
 1001.1637  2041.1368  1716.1817  1485.7563  2379.8347  2105.8293  1788.2276 
sd.a.gl[8] sd.a.gl[9] 
 2878.3639   828.0619 
Show code
node.sub <- grep("^R\\[", colnames(silv.samples[[1]]), value = TRUE)
silv.samples[, node.sub[sample(1:length(node.sub), 36)], drop = FALSE] %>%
mcmc_trace()

Show code
pl <- silv.samples[, node.sub[], drop = FALSE] %>%
gelman.diag(multivariate = FALSE) # Error in chol.default(W) : the leading minor of order 1 is not positive
pl <- silv.samples[, node.sub[], drop = FALSE] %>%
effectiveSize()
print("effective sample size")
[1] "effective sample size"
Show code
pl[which(pl < 1000 & pl > 0 )]
 R[1, 3]  R[1, 4]  R[3, 4] 
799.4678 790.3968 984.1632 

Silvering variables

Show code
silv.samples[[1]] %>%
  gather_draws(mu.a.gl[er], sep = ",") %>%
  ungroup() %>%
  left_join(data.silv %>% distinct(er,ecoreg)) %>%
  #filter(!is.na(ecoreg)) %>%
  ggplot() +
  geom_density(aes(x = exp(.value), color = ecoreg)) +
  facet_wrap(~ecoreg, scales = "free_y", ncol = 3) +
  theme_minimal() +
  labs(title = "mu.a.gl")
ungroup: no grouping variables remain
distinct: removed 147,485 rows (>99%), 9 rows remaining
Joining with `by = join_by(er)`
left_join: added one column (ecoreg)
           > rows only in x                               0
           > rows only in data.silv %>% distinct(.. (     0)
           > matched rows                            13,500
           >                                        ========
           > rows total                              13,500

Show code
silv.samples[[1]] %>%
  gather_draws(alpha[er, yr], sep = ",") %>%
  ungroup() %>%
  left_join(data.silv %>% distinct(er,ecoreg)) %>%
  ggplot() +
  geom_density(aes(x = .value, color = factor(yr))) +
  facet_wrap(~ecoreg, scales = "free_y", ncol = 3) +
  theme_minimal() +
  labs(title = "alpha")
ungroup: no grouping variables remain
distinct: removed 147,485 rows (>99%), 9 rows remaining
Joining with `by = join_by(er)`
left_join: added one column (ecoreg)
           > rows only in x                                0
           > rows only in data.silv %>% distinct(.. (      0)
           > matched rows                            378,000
           >                                        =========
           > rows total                              378,000

Show code
silv.samples[[1]] %>%
  gather_draws(bs[er], sep = ",") %>%
  ungroup() %>%
  left_join(data.silv %>% distinct(er,ecoreg)) %>%
  ggplot() +
  geom_density(aes(x = .value, color = ecoreg)) +
  theme_minimal() +
  labs(title = "bs")
ungroup: no grouping variables remain
distinct: removed 147,485 rows (>99%), 9 rows remaining
Joining with `by = join_by(er)`
left_join: added one column (ecoreg)
           > rows only in x                               0
           > rows only in data.silv %>% distinct(.. (     0)
           > matched rows                            13,500
           >                                        ========
           > rows total                              13,500

Show code
silv.samples[[1]] %>%
  gather_draws(s.sig,bt) %>%
  ggplot() +
  geom_density(aes(x = .value, color = .variable)) +
  facet_wrap(~.variable, scales = "free") +
  theme_minimal() 

Predicted / observed

Show code
dfp <- silv.samples[[1]] %>%
  spread_draws(s.mu[i]) %>%
  median_qi() %>%
  rename(pred.sl = s.mu,
         ind.id = i) %>%
  ungroup() %>%
  left_join(data.silv %>% 
              select(length,ecoreg) %>% 
              rename(obs.sl = length) %>%
              mutate(ind.id = row_number()))
rename: renamed 2 variables (ind.id, pred.sl)
ungroup: no grouping variables remain
select: dropped 349 variables (mei_fi_id, fi_year, fi_date, ser_x, ser_y, …)
rename: renamed one variable (obs.sl)
mutate: new variable 'ind.id' (integer) with 147,494 unique values and 0% NA
Joining with `by = join_by(ind.id)`
left_join: added 2 columns (obs.sl, ecoreg)
           > rows only in x                                0
           > rows only in data.silv %>% select(le.. (      0)
           > matched rows                            147,494
           >                                        =========
           > rows total                              147,494
Show code
dfp %>%
  slice_sample(n = 5000) %>%
  ggplot() +
  geom_abline(slope =1, color = "darkgrey") +
  geom_point(aes(x = obs.sl, y = exp(pred.sl)), shape = 4,alpha = 0.7) +
  facet_wrap(~ecoreg) +
  theme_minimal() 
slice_sample: removed 142,494 rows (97%), 5,000 rows remaining

Show code
ggsave(paste0(home, "/report/silv_po.png"), width = 17, height = 14, units = "cm")

Length over time

Show code
dfp <- silv.samples[[1]] %>%
  spread_draws(alpha[er, yr], bs[er],sep = ",") %>%
  mutate(male = alpha + bs,
         year = yr + 1997) %>%
  rename(female = alpha) %>%
  pivot_longer(cols = c("female","male"), names_to = "sex", values_to = "length") %>%
  mutate(length = exp(length)) %>%
  left_join(data.silv %>% distinct(er,ecoreg)) %>%
  ungroup() %>%
  mutate(m_length = median(length), .by = c(year,sex))
mutate (grouped): new variable 'male' (double) with 377,496 unique values and 0% NA
                  new variable 'year' (double) with 28 unique values and 0% NA
rename: renamed one variable (female)
pivot_longer: reorganized (female, male) into (sex, length) [was 378000x9, now 756000x9]
mutate (grouped): changed 756,000 values (100%) of 'length' (0 new NAs)
distinct: removed 147,485 rows (>99%), 9 rows remaining
Joining with `by = join_by(er)`
left_join: added one column (ecoreg)
           > rows only in x                                0
           > rows only in data.silv %>% distinct(.. (      0)
           > matched rows                            756,000
           >                                        =========
           > rows total                              756,000
ungroup: no grouping variables remain
mutate: new variable 'm_length' (double) with 56 unique values and 0% NA
Show code
dfp %>%
  filter(sex == "female") %>%
  ungroup() %>%
  ggplot() +
  geom_boxplot(aes(y = length, x = year, group = year), color = "#F8766D", outliers = FALSE) +
  geom_line(aes(y = m_length, x = year), color = "#F8766D", alpha = 0.4) +
  facet_wrap(~ecoreg, scales = "free_y") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = rel(.75))) +
  labs(x = "year")
filter: removed 378,000 rows (50%), 378,000 rows remaining
ungroup: no grouping variables remain

Show code
ggsave(paste0(home, "/report/silv_lf.png"), width = 17, height = 14, units = "cm")
  
dfp %>%
  filter(sex == "male") %>%
  ungroup() %>%
  ggplot() +
  geom_boxplot(aes(y = length, x = year, group = year), color = "#00BFC4", outliers = FALSE) +
  geom_line(aes(y = m_length, x = year), color = "#00BFC4", alpha = 0.4) +
  scale_x_discrete(breaks = function(x) x[as.numeric(x) %% 2 == 0]) +
  facet_wrap(~ecoreg, scales = "free_y") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, size = rel(.75))) +
  labs(x = "year")
filter: removed 378,000 rows (50%), 378,000 rows remaining
ungroup: no grouping variables remain

Show code
ggsave(paste0(home, "/report/silv_lm.png"), width = 17, height = 14, units = "cm")

silv.samples %>%
  spread_draws(alpha[er, yr], bs[er],sep = ",") %>%
  median_qi() %>%
  mutate(male = alpha + bs,
         year = yr + 1997) %>%
  rename(female = alpha) %>%
  pivot_longer(cols = c("female","male"), names_to = "sex", values_to = "length") %>%
  left_join(data.silv %>% distinct(er,ecoreg)) %>%
  mutate(length_mc = exp(length) - mean(exp(length)), .by = c(sex)) %>%
  ggplot() +
  geom_line(aes(y = length_mc, x = year, color = ecoreg)) +
  facet_wrap(~sex, scales = "free_y") +
  theme_minimal() +
  labs(x = "year", caption = "centered by mean by sex")
mutate: new variable 'male' (double) with 252 unique values and 0% NA
        new variable 'year' (double) with 28 unique values and 0% NA
rename: renamed one variable (female)
pivot_longer: reorganized (female, male) into (sex, length) [was 252x13, now 504x13]
distinct: removed 147,485 rows (>99%), 9 rows remaining
Joining with `by = join_by(er)`
left_join: added one column (ecoreg)
           > rows only in x                            0
           > rows only in data.silv %>% distinct(.. (  0)
           > matched rows                            504
           >                                        =====
           > rows total                              504
mutate: new variable 'length_mc' (double) with 504 unique values and 0% NA

Correlation/covariance matrix

Show code
# median sample correlation 
samp <- silv.samples[[1]]

cols <- grep("^R",colnames(as.matrix(samp)))


corr_arr <- array(NA, dim = c(9, 9, nrow(samp)))

for (i in 1:nrow(samp)) {

  R_sample <- matrix(
    as.numeric(samp[i, cols]),
    nrow = 9,
    ncol = 9,
    byrow = FALSE
  )

  corr_arr[, , i] <- crossprod(R_sample)
}

corr <- apply(corr_arr, c(1, 2), median)

var.lats <- data.silv %>%
  mutate(m.lat = mean(ser_y), .by = ecoreg) %>%
  distinct(er, ecoreg, m.lat) %>%
  arrange(-m.lat)
mutate: new variable 'm.lat' (double) with 9 unique values and 0% NA
distinct: removed 147,485 rows (>99%), 9 rows remaining
Show code
er.labs <- var.lats %>%
  select(er, ecoreg) %>%
  { setNames(.$ecoreg, .$er) }
select: dropped one variable (m.lat)
Show code
as.data.frame(corr) %>%
  mutate(r = row_number()) %>%
  pivot_longer(cols = -r, names_to = "c", values_to = "value") %>%
  mutate(c = as.integer(gsub("V", "", c)) ) %>%
  left_join(var.lats, by = c("r" = "er")) %>%
  left_join(var.lats, by = c("c" = "er"), suffix = c(".x", ".y")) %>%
  mutate(r = reorder(r, m.lat.x),
         c = reorder(c, -m.lat.y)) %>%
  ggplot(aes(x=c, y=r, fill = value)) + 
  geom_tile() +
  scale_fill_viridis_c() +
  scale_x_discrete(labels = er.labs) +
  scale_y_discrete(labels = er.labs) +
  theme(axis.text = element_text(angle = 45, hjust = 1, size = rel(.75)))
mutate: new variable 'r' (integer) with 9 unique values and 0% NA
pivot_longer: reorganized (V1, V2, V3, V4, V5, …) into (c, value) [was 9x10, now 81x3]
mutate: converted 'c' from character to integer (0 new NA)
left_join: added 2 columns (ecoreg, m.lat)
           > rows only in x          0
           > rows only in var.lats ( 0)
           > matched rows           81
           >                       ====
           > rows total             81
left_join: added 4 columns (ecoreg.x, m.lat.x, ecoreg.y, m.lat.y)
           > rows only in x          0
           > rows only in var.lats ( 0)
           > matched rows           81
           >                       ====
           > rows total             81
mutate: converted 'r' from integer to factor (0 new NA)
        converted 'c' from integer to factor (0 new NA)

Show code
# var.temps <- data.silv %>%
#   mutate(m.temp = mean(temp.sc), .by = ecoreg) %>%
#   distinct(er, ecoreg, m.temp) %>%
#   arrange(-m.temp)
# 
# as.data.frame(corr) %>%
#   mutate(r = row_number()) %>%
#   pivot_longer(cols = -r, names_to = "c", values_to = "value") %>%
#   mutate(c = as.integer(gsub("V", "", c)) ) %>%
#   left_join(var.temps, by = c("r" = "er")) %>%
#   left_join(var.temps, by = c("c" = "er"), suffix = c(".x", ".y")) %>%
#   mutate(r = reorder(r, m.temp.x),
#          c = reorder(c, -m.temp.y)) %>%
#   ggplot(aes(x=c, y=r, fill = value)) + 
#   geom_tile() +
#   scale_fill_viridis_c() +
#   scale_x_discrete(labels = er.labs) +
#   scale_y_discrete(labels = er.labs) +
#   theme(axis.text = element_text(angle = 45, hjust = 1, size = rel(.75)))

# one sample alternative
# R_sample <- matrix(samp[nrow(samp), cols], 9, 9)
# corr = crossprod(R_sample)  ## i.e., t(R_sample) %*% R_sample
# melt(corr) %>%
#   left_join(var.lats, by = c("Var1" = "er")) %>%
#   left_join(var.lats, by = c("Var2" = "er"), suffix = c(".x", ".y")) %>%
#   mutate(Var1 = reorder(Var1, m.lat.x),
#          Var2 = reorder(Var2, -m.lat.y)) %>%
#   ggplot(aes(x=Var2, y=Var1, fill = value)) + 
#   geom_tile() +
#   scale_fill_viridis_c() +
#   scale_x_discrete(labels = er.labs) +
#   scale_y_discrete(labels = er.labs) +
#   theme(axis.text = element_text(angle = 45, hjust = 1, size = rel(.75)))